#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include <string.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sort_long.h>
#ifdef USEOPENMP
#include <omp.h>
#endif

#include "halofinder.h"

int main(int argc, char **argv)
{
  /* vars */
  char name[MAX_FILENAME];
#ifdef USEOPENMP
  int provided;
#endif
  
  /* init MPI and get current tasks and number of tasks */
  /*#ifdef USEOPENMP
  int rc = MPI_Init_thread(&argc,&argv,MPI_THREAD_FUNNELED,&provided);
  if(rc != MPI_SUCCESS || provided != MPI_THREAD_FUNNELED)
  {
  fprintf(stderr,"Error starting MPI program. Terminating.\n");
  MPI_Abort(MPI_COMM_WORLD,rc);
  }    
  #else*/
  int rc = MPI_Init(&argc,&argv);
  if(rc != MPI_SUCCESS)  
    {
      fprintf(stderr,"Error starting MPI program. Terminating.\n");
      MPI_Abort(MPI_COMM_WORLD,rc);
    }
  //#endif
  
  MPI_Comm_size(MPI_COMM_WORLD,&NTasks);
  MPI_Comm_rank(MPI_COMM_WORLD,&ThisTask);
  
  logProfileTag(PROFILETAG_TOTTIME);
  
#ifdef MEMWATCH
  mwDoFlush(1);
#endif
  
  /* init all vars for raytracing 
     1) read from config file
     2) set search radius
     3) set plane names
  */
  read_config_file(argv[1]); fflush(stderr);
  
  //test for proper number of threads
#ifdef USEOPENMP
  omp_set_num_threads(atoi(argv[2]));
  if(ThisTask == 0)
    {
      fprintf(stderr,"%d: # of tasks = %d, # of threads = %d, requested %d threads\n",ThisTask,NTasks,omp_get_max_threads(),atoi(argv[2]));
      fflush(stderr);
#pragma omp parallel
      {
	fprintf(stderr,"%d: thread %d\n",ThisTask,omp_get_thread_num());
	fflush(stderr);
      }
    }
  fflush(stderr);
#endif
  
#ifdef TEST_CODE
  test_shiftpartshalos();
  //test_kdtree_nbrsfind();
  //test_kdtree_knnbrsfind();
  test_fofhalogroupfinder();
  //test_sort_halos();
  //test_peakhaloio();
  //test_scattergatherpeakshalos();
  

  //test_exchangebufferregions();
#else
  logProfileTag(PROFILETAG_DOMAINDECOMP);
  setupdomaindecomp();
  logProfileTag(PROFILETAG_DOMAINDECOMP);
  
  //read particles
  logProfileTag(PROFILETAG_PARTIO);
  readpartsfromsnapshot();
  logProfileTag(PROFILETAG_PARTIO);
  
  logProfileTag(PROFILETAG_DOMAINDECOMP);
  if(NumDomains[0] > 1)
    exchangebufferregions();
  if(NumDomains[0] > 1)                /* shift all parts and halos to domain coords - won't shift back unitl later */
    shiftpartshalosglobal2domain();
  reorder_parts_tree();
  build_parts_kdtree();                /* building main kdtree for nearest neighbors searching */
  logProfileTag(PROFILETAG_DOMAINDECOMP);
  
  //find peaks
  if(strlen(haloFinderData.PeakInputFile) == 0)
    {
      logProfileTag(PROFILETAG_PEAKFIND);
      peakfinder();
      logProfileTag(PROFILETAG_PEAKFIND);
      fprintf(stderr,"%d: peak find time = %lf [s]\n",ThisTask,getTotTimeProfileTag(PROFILETAG_PEAKFIND));
      
#ifdef PERCOLATEPEAKS
      logProfileTag(PROFILETAG_PEAKPERC);
      peakpercolation();
      logProfileTag(PROFILETAG_PEAKPERC);
#endif
    }
  else
    {
      logProfileTag(PROFILETAG_PEAKHALOIO);
      readpeaks();
      if(NumDomains[0] > 1)                /* after reading pekas need ot shift coords to domain */
	shifthalosglobal2domain();
      logProfileTag(PROFILETAG_PEAKHALOIO);
    }
  
  ///////////////////////////
  // make peaks into halos //
  ///////////////////////////
  logProfileTag(PROFILETAG_HALORAD);
  computehaloradii();
  logProfileTag(PROFILETAG_HALORAD);
  fprintf(stderr,"%d: halo radius compute time = %lf [s]\n",ThisTask,getTotTimeProfileTag(PROFILETAG_HALORAD));
  
  //write peaks
  if(strlen(haloFinderData.PeakOutputFile) != 0)
    {
      logProfileTag(PROFILETAG_PEAKHALOIO);
      if(NumDomains[0] > 1)                
	shifthalosdomain2global();
      writepeaks();
      if(NumDomains[0] > 1)                
	shifthalosglobal2domain();
      logProfileTag(PROFILETAG_PEAKHALOIO);
    }
  
  logProfileTag(PROFILETAG_DOMAINDECOMP);
  destroyparts();
  destroy_parts_kdtree();
  logProfileTag(PROFILETAG_DOMAINDECOMP);
  
  //timing tags are in function calls
#ifndef PARALLELHALOPERC
  if(NumDomains[0] > 1)
    shifthalosdomain2global();
  gatherpeakshalos();
  halopercolation_serial();
#else
  halopercolation_parallel();
  if(NumDomains[0] > 1)
    shifthalosdomain2global();
#endif
  
  //used to get timing infor correct for halo i/o call
  MPI_Barrier(MPI_COMM_WORLD);
  
  logProfileTag(PROFILETAG_PEAKHALOIO);
  writehalos();
  logProfileTag(PROFILETAG_PEAKHALOIO);
  
  logProfileTag(PROFILETAG_DOMAINDECOMP);
  destroyhalos();
  destroyparts();
  destroydomaindecomp();
  logProfileTag(PROFILETAG_DOMAINDECOMP);
#endif /* TEST_CODE */  
  
  /* finish profiling info*/
  logProfileTag(PROFILETAG_TOTTIME);
  sprintf(name,"./timing");
  const char *ProfileTagNames[] = {"totime","partio","peakfind","peakhaloio","halorad","haloperc","peakperc","domaindecomp"};
  printProfileInfo(name,ProfileTagNames);
  
  MPI_Finalize();
  return 0;
}
